Lets start by loading in libraries we will need
library("xlsx")
library(ggplot2)
library(forcats)
library(tidyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(priceTools)
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
library(dummies)
## dummies-1.5.6 provided by Decision Patterns
and the data we will be using.
isotope_data <- read.xlsx("edited_community_isotope_data.xlsx",2)
Simple visualization of the data
head(isotope_data)
## Panel_Site Panel_Number Site_and_Panel Species_ID Taxa
## 1 FC 1 FC_01 Alcyonidium polypylum Bryozoan
## 2 FC 1 FC_01 Alcyonidium polypylum Bryozoan
## 3 FC 1 FC_01 Alitta succinea worm
## 4 FC 1 FC_01 Amphibalanus eburneus Barnacle
## 5 FC 1 FC_01 Bugula neritina Bryozoan
## 6 FC 1 FC_01 Conopeum chesapeakensis Bryozoan
## pedino_13C_enrichment atom_13C_percent pico_15N_enrichment atom_15N_percent
## 1 -18.113866 1.085858 6.193911 0.3685638
## 2 -20.589898 1.083149 6.303070 0.3686036
## 3 -21.085582 1.082607 12.824082 0.3709834
## 4 -19.677066 1.084148 7.408607 0.3690071
## 5 2.423247 1.108318 6.323871 0.3686112
## 6 -7.060286 1.097948 8.278584 0.3693246
unique(isotope_data$Taxa)
## [1] "Bryozoan" "worm" "Barnacle" "Bivalve" "Amphipod"
## [6] "Anemone" "Tunicate" "brittle star" "sponge" "crab"
## [11] "Isopod"
The values contained in the columns pedino_13C_enrichment and pico_15N_enrichment represent FILL IN
The values contained in the columns atom_13C_percent and atom_15N_percent represent FILL IN
Barplots of the Means across each column of interest for Taxa
#barplots by taxa means
par(mfrow=c(2,2))
color_vector_taxa = rev(rainbow(length(unique(isotope_data$Taxa))))
mean_taxa_pedino = tapply(isotope_data$pedino_13C_enrichment, isotope_data$Taxa, mean)
mean_taxa_pedino_percent = tapply(isotope_data$atom_13C_percent,isotope_data$Taxa, mean)
barplot(mean_taxa_pedino,las=2, col = color_vector_taxa, ylab = 'Pedino Enrichment (13C) Mean', main = "Pedino Enrichment Mean vs. Taxa")
barplot(mean_taxa_pedino_percent, las=2, ylim=c(0.0, 1.2), col = color_vector_taxa, ylab = 'Pedino 13C atom% mean', main = "Pedino Atom% vs. Taxa")
mean_taxa_pico = tapply(isotope_data$pico_15N_enrichment, isotope_data$Taxa, mean)
mean_taxa_pico_percent = tapply(isotope_data$atom_15N_percent,isotope_data$Taxa, mean)
barplot(mean_taxa_pico,las=2, col = color_vector_taxa, ylab = "Pico Enrichment (15N) Mean", main = "Pico Enrichment Mean vs. Taxa")
barplot(mean_taxa_pico_percent, las=2, col = color_vector_taxa,ylab = "Pico Atom% (15N)", main = "Pico Atom% Mean vs. Taxa")
Barplots of the Means across variable each column for Site
#barplots by site
par(mfrow=c(2,2))
color_vector_site = rev(terrain.colors(length(unique(isotope_data$Taxa))))
mean_site_pedino = tapply(isotope_data$pedino_13C_enrichment, isotope_data$Panel_Site, mean)
mean_site_pedino_percent = tapply(isotope_data$atom_13C_percent,isotope_data$Panel_Site, mean)
barplot(mean_site_pedino,las=2, col = color_vector_site, ylab = "Pedino Enrichment (13C) Mean", main = "Pedino Enrichment vs Site")
barplot(mean_site_pedino_percent, las=2, ylim=c(0.0, 1.2), col = color_vector_site,ylab = "Pedino Atom% (13C)", main = "Pedino Atom% Mean vs. Site")
mean_site_pico = tapply(isotope_data$pico_15N_enrichment, isotope_data$Panel_Site, mean)
mean_site_pico_percent = tapply(isotope_data$atom_15N_percent,isotope_data$Panel_Site, mean)
barplot(mean_site_pico,las=2, col = color_vector_site, ylab = "Pico Enrichment (15N) Mean", main = "Pico Enrichment Mean vs. Site")
barplot(mean_site_pico_percent, las=2, col = color_vector_site,ylab = "Pico Atom% (15N) Mean", main = "Pico Atom% Mean vs. Site")
Now we will make some boxplots for Taxa and Enrichment
ggplot(data = isotope_data) +
geom_boxplot(mapping = aes(x = Taxa, y = pedino_13C_enrichment, color = Taxa)) +
labs(x = 'Taxa', y = 'pedino_13C_enrichment', title = 'taxa vs pedino_13C_enrichment') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = isotope_data) +
geom_boxplot(mapping = aes(x = Taxa, y = pico_15N_enrichment, color = Taxa)) + labs(x = 'Taxa', y = 'pico_15N_enrichment', title = 'taxa vs pico_15N_enrichment') + theme(axis.text.x = element_text(angle = 90))
As well as boxplots for Taxa and Atom %
#Pedino
ggplot(data = isotope_data) +
geom_boxplot(mapping = aes(x = Taxa, y = atom_13C_percent, color = Taxa)) +
labs(x = 'Taxa', y = 'atom 13C %', title = 'taxa vs pedino 13C %') + theme(axis.text.x = element_text(angle = 90))
#pico
ggplot(data = isotope_data) +
geom_boxplot(mapping = aes(x = Taxa, y = atom_15N_percent, color = Taxa)) +
labs(x = 'Taxa', y = "atom 15N %", title = 'taxa vs pico 15N %') + theme(axis.text.x = element_text(angle = 90))
basic taxa linear regression pedino
pedino_taxa_lm = lm(isotope_data$pedino_13C_enrichment ~ -1 + isotope_data$Taxa)
summary(pedino_taxa_lm)
##
## Call:
## lm(formula = isotope_data$pedino_13C_enrichment ~ -1 + isotope_data$Taxa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.363 -2.769 -0.437 1.206 182.301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## isotope_data$TaxaAmphipod 8.132 3.538 2.298 0.0222 *
## isotope_data$TaxaAnemone -18.618 7.222 -2.578 0.0104 *
## isotope_data$TaxaBarnacle -20.353 1.479 -13.757 < 2e-16 ***
## isotope_data$TaxaBivalve 29.146 5.594 5.210 3.46e-07 ***
## isotope_data$Taxabrittle star -9.580 12.509 -0.766 0.4444
## isotope_data$TaxaBryozoan -3.910 2.553 -1.531 0.1267
## isotope_data$Taxacrab -6.921 12.509 -0.553 0.5805
## isotope_data$TaxaIsopod -1.754 5.107 -0.343 0.7315
## isotope_data$Taxasponge 36.403 8.845 4.115 4.96e-05 ***
## isotope_data$TaxaTunicate -3.953 3.127 -1.264 0.2072
## isotope_data$Taxaworm -5.502 2.948 -1.866 0.0630 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.69 on 309 degrees of freedom
## Multiple R-squared: 0.4509, Adjusted R-squared: 0.4313
## F-statistic: 23.06 on 11 and 309 DF, p-value: < 2.2e-16
basic taxa linear regression pico
pico_taxa_lm = lm(isotope_data$pico_15N_enrichment ~ -1 + isotope_data$Taxa)
summary(pico_taxa_lm)
##
## Call:
## lm(formula = isotope_data$pico_15N_enrichment ~ -1 + isotope_data$Taxa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3637 -1.6963 -0.2294 0.8359 24.0512
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## isotope_data$TaxaAmphipod 9.5808 0.5611 17.075 < 2e-16 ***
## isotope_data$TaxaAnemone 10.7529 1.1453 9.388 < 2e-16 ***
## isotope_data$TaxaBarnacle 9.5564 0.2346 40.734 < 2e-16 ***
## isotope_data$TaxaBivalve 10.4804 0.8872 11.813 < 2e-16 ***
## isotope_data$Taxabrittle star 8.3772 1.9838 4.223 3.18e-05 ***
## isotope_data$TaxaBryozoan 8.7510 0.4049 21.611 < 2e-16 ***
## isotope_data$Taxacrab 10.0385 1.9838 5.060 7.20e-07 ***
## isotope_data$TaxaIsopod 5.6950 0.8099 7.032 1.31e-11 ***
## isotope_data$Taxasponge 12.6813 1.4027 9.040 < 2e-16 ***
## isotope_data$TaxaTunicate 9.5369 0.4959 19.230 < 2e-16 ***
## isotope_data$Taxaworm 11.1417 0.4676 23.828 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.805 on 309 degrees of freedom
## Multiple R-squared: 0.924, Adjusted R-squared: 0.9213
## F-statistic: 341.6 on 11 and 309 DF, p-value: < 2.2e-16
Boxplots for enrichment by taxa mean- per panel
FC
FC_isotope <- subset(isotope_data, isotope_data$Panel_Site == "FC")
ggplot(data = FC_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment FC') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = FC_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment FC') + theme(axis.text.x = element_text(angle = 90))
HI
HI_isotope <- subset(isotope_data, isotope_data$Panel_Site == "HI")
ggplot(data = HI_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment HI') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = HI_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment HI') + theme(axis.text.x = element_text(angle = 90))
ID
ID_isotope <- subset(isotope_data, isotope_data$Panel_Site == "ID")
ggplot(data = ID_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment ID') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = ID_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment ID') + theme(axis.text.x = element_text(angle = 90))
IRL_3
IRL3_isotope <- subset(isotope_data, isotope_data$Panel_Site == "IRL3")
ggplot(data = IRL3_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment IRL3') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = IRL3_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment IRL3') + theme(axis.text.x = element_text(angle = 90))
MIM
MIM_isotope <- subset(isotope_data, isotope_data$Panel_Site == "MIM")
ggplot(data = MIM_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment MIM') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = MIM_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment MIM') + theme(axis.text.x = element_text(angle = 90))
MO2
MO2_isotope <- subset(isotope_data, isotope_data$Panel_Site == "MO2")
ggplot(data = MO2_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment MO2') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = MO2_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment MO2') + theme(axis.text.x = element_text(angle = 90))
SCD
SCD_isotope <- subset(isotope_data, isotope_data$Panel_Site == "SCD")
ggplot(data = SCD_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment SCD') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = SCD_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment SCD') + theme(axis.text.x = element_text(angle = 90))
SMS
SMS_isotope <- subset(isotope_data, isotope_data$Panel_Site == "SMS")
ggplot(data = SMS_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment SMS') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = SMS_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment SMS') + theme(axis.text.x = element_text(angle = 90))
WP
WP_isotope <- subset(isotope_data, isotope_data$Panel_Site == "WP")
ggplot(data = WP_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment WP') + theme(axis.text.x = element_text(angle = 90))
ggplot(data = WP_isotope) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment WP') + theme(axis.text.x = element_text(angle = 90))
Now we will create a dummy matrix (isotope_copy) in order to reorder the factor levels of Species ID by Taxa to graph a master plot of the Species and Enrichment values in an organized manner
isotope_copy <- isotope_data
isotope_copy$Species_ID <- factor(isotope_copy$Species_ID)
isotope_copy$Species_ID <- fct_reorder(isotope_copy$Species_ID, isotope_copy$Taxa, .fun = max)
Species versus Pedino Enrichment
ggplot(data = isotope_copy) +
geom_boxplot(mapping = aes(x = Species_ID, y = pedino_13C_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pedino_13C_enrichment', title = 'species vs pedino_13C_enrichment master') + theme(axis.text.x = element_text(angle = 90))
Species versus Pico enrichment
ggplot(data = isotope_copy) +
geom_boxplot(mapping = aes(x = Species_ID, y = pico_15N_enrichment, color=Taxa)) +
labs(x = 'species_id', y = 'pico_15N_enrichment', title = 'species vs pico_15N_enrichment master') + theme(axis.text.x = element_text(angle = 90))
subsetting the species of algae
pedino_subset <- isotope_data[c('Site_and_Panel',"Species_ID","pedino_13C_enrichment","atom_13C_percent")]
pico_subset <- isotope_data[c('Site_and_Panel',"Species_ID","pico_15N_enrichment", "atom_15N_percent")]
creating community matrixes
comm_pedino <- with(pedino_subset, tapply(pedino_13C_enrichment, list(Site_and_Panel, Species_ID), sum))
comm_pedino <- ifelse(is.na(comm_pedino), 0, comm_pedino)
comm_pico <- with(pico_subset, tapply(pico_15N_enrichment, list(Site_and_Panel, Species_ID), sum))
comm_pico <- ifelse(is.na(comm_pico), 0, comm_pico)
pico pca
pico_pca = rda(comm_pico, scale = TRUE)
pico_pca
## Call: rda(X = comm_pico, scale = TRUE)
##
## Inertia Rank
## Total 46
## Unconstrained 46 25
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 7.845 6.618 4.749 3.905 3.617 2.546 2.308 2.109
## (Showing 8 of 25 unconstrained eigenvalues)
# the eigenvalues sum up to equal the total interia (i.e., total variance in this case)
sum(pico_pca$CA$eig)
## [1] 46
# the ratio of the eigenvalue to the total variance is the amount of
# variance explained by each PCA axis
round(pico_pca$CA$eig / pico_pca$tot.chi, 2)
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16
## 0.17 0.14 0.10 0.08 0.08 0.06 0.05 0.05 0.04 0.04 0.03 0.03 0.02 0.02 0.02 0.02
## PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25
## 0.01 0.01 0.01 0.01 0.01 0.00 0.00 0.00 0.00
#We can see from above that the PCA axis 1 captures approximately 17% of the total variance in the community matrix. Let's graph the data to better get a sense of the correlation structure.
plot(pico_pca)
biplot(pico_pca)
ordiplot(pico_pca, display = 'species')
orditorp(pico_pca, display = 'species')
pedino pca
pedino_pca = rda(comm_pedino, scale = TRUE)
pedino_pca
## Call: rda(X = comm_pedino, scale = TRUE)
##
## Inertia Rank
## Total 46
## Unconstrained 46 25
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 7.324 6.879 5.468 3.627 3.328 2.813 2.466 2.161
## (Showing 8 of 25 unconstrained eigenvalues)
par(mfrow= c(1,1))
# the eigenvalues sum up to equal the total interia (i.e., total variance in this case)
sum(pedino_pca$CA$eig)
## [1] 46
# the ratio of the eigenvalue to the total variance is the amount of
# variance explained by each PCA axis
round(pedino_pca$CA$eig / pedino_pca$tot.chi, 2)
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16
## 0.16 0.15 0.12 0.08 0.07 0.06 0.05 0.05 0.04 0.03 0.03 0.03 0.02 0.02 0.02 0.02
## PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25
## 0.01 0.01 0.01 0.01 0.00 0.00 0.00 0.00 0.00
#We can see from above that the PCA axis 1 captures approximately 16% of the total variance in the community matrix. Let's graph the data to better get a sense of the correlation structure.
plot(pedino_pca)
biplot(pedino_pca)
ordiplot(pedino_pca, display = 'species')
orditorp(pedino_pca, display = 'species')